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We present a Quantum Monte Carlo study of the dissociation energy and the dispersion curve of the water 
dimer, a prototype of hydrogen bonded system. Our calculations are based on a wave function which is a 
modern and fully correlated implementation of the Pauling's valence bond idea: the Jastrow Antisymmetrised 
Geminal Power (JAGP) [Casula et al J.Chem.Phys. 1 19,6500-65 1 1 (2003)]. With this variational wave function 
we obtain a binding energy of —4.5(0.1) kcal/mol, that is only slightly increased to —4.9(0.1) kcal/mol by 
using the Lattice Regularized Diffusion Monte Carlo (LRDMC). This projection technique allows to improve 
substantially the correlation energy of a given variational guess and indeed, when applied to the JAGP, yields a 
binding energy in fair agreement with the value of —5.0 kcal/mol reported by experiments and other theoretical 
works. The minimum position, the curvature and the asymptotic behavior of the dispersion curve are well 
reproduced both at the variational and LRDMC level. Moreover, thanks to the simplicity and the accuracy of 
our variational approach, we are able to dissect the various contributions to the binding energy of the water dimer 
in a systematic and controlled way. This is achieved by appropriately switching off determinantal and Jastrow 
variational terms in the JAGP. Within this scheme, we estimate that the van der Waals contribution to the electron 
correlation is substantial and amounts to 1.5(0.2) kcal/mol, this value being comparable with the intermolecular 
covalent energy, that we find to be 1.1(0.2) kcal/mol. The present Quantum Monte Carlo approach based on the 
JAGP wavefunction reveals as a promising tool for the interpretation and the quantitative description of weakly 
interacting systems, where both dispersive and covalent energy contributions play an important role. 



PACS numbers: 

I. INTRODUCTION 

The hydrogen bond is a fundamental intramolecular and in- 
termolecular interaction determining the properties of a large 
number of systems from liquids to solids, from biological [ 1 ] 
to inorganic yfl. Hydrogen bond is commonly defined as a 
local bond in which an hydrogen atom is attached to an elec- 
tronegative group (the donor) interacting with another nearby 
electronegative group (the acceptor), not covalently attached. 
Dissociation energies cover a range of about two orders of 
magnitude, ranging from —0.2 to 40 kcal/mol, the H-bonding 
arising from the interplay of different types of interactions. 
Electrostatic forces play the major role in a large number of 
hydrogen bonds, although charge transfer effects and van der 
Waals (vdW) interactions are always present. 

Water, the most studied H-bonding liquid, represents the 
prototype of hydrogen bonded systems. The energetics and 
directionality of water hydrogen bond is a key factor for un- 
derstanding the anomalous properties of the water phase di- 
agram [3], the behavior of small water clusters 0. El 0] and 
the role of aqueous environment in a variety of biological sys- 
tems \M. The dissociation energy of the isolated water dimer 
lies in the middle of the hydrogen bond dissociation energy 
scale, the most stable configuration being associated with a 
binding energy £)fff = —5.0 kcal/mol, as extrapolated by 
experimental data [8]. The partitioning of this energy in dif- 
ferent contribution terms is still subject of vivid debate. At 
the equilibrium bonding distance, typically in a range between 



2.5 and 3.5 A, quantum effects become relevants and a pure 
electrostatic picture of the interaction is not fully satisfactory. 
The partial covalent nature of the hydrogen bond has been re- 
cently invoked by a first analysis of the Compton profile on 
ice Ih 091]. However, the interpretation of the experimental 
data has been questioned by several authors ITol IllH and fur- 
ther revised 1121 . The amount of the intermolecular covalent 
contribution, if any, to the binding energy is still an open is- 
sue. On the other hand, due to the lack of an unambiguous 
computational protocol it is still not clear how to estimate the 
van der Waals contribution to the hydrogen bonding. The role 
of these interactions may also be at the basis of the current 
drawbacks of empirical forcefields in use for large scale sim- 
ulations iflUflUl. 

A definition of the intermolecular covalent component of 
the hydrogen bonding can be drawn using the intuitive pic- 
ture of chemical bond introduced by Pauling as the superposi- 
tion of Lewis' structures [15]. In the simple case of hydrogen 
bonding in a water dimer (^0)2, three mesomeric Lewis 
structures may be drawn, one of them describing the charge 
transfer situation {OH)~ ■ ■ ■ (H^O) + ', that confers partial co- 
valent character to the hydrogen bond. Within a quantitative 
Valence Bond representation it would be therefore possible to 
distinguish, in a simple way, the covalent intermolecular en- 
ergy contribution by the other interaction energy terms. 

At the same time, because of the crucial role of electronic 
correlation, especially for dispersive interactions, high quality 
electronic structure correlated methods (based on molecular 
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orbital theory) are necessary to a proper quantitative descrip- 
tion of hydrogen bond. New classes of Density Functional 
Theory (DFT) functionals have been developed in the past 
years with the aim to describe weakly bound systems avoid- 
ing semi empirical approaches. Nevertheless, the highly non 
perturbative and non-local character of the vdW interactions 
makes difficult their inclusions in DFT schemes without the 
resorting to ad hoc empirical parameterizations II 1611 . Look- 
ing for an ab initio method free of empiricism, the Quan- 
tum Monte Carlo appears as a possible alternative 
to other more standard quantum chemistry methods such as 
Configuration Interaction, Moller Plesset perturbation theory 
or coupled-cluster (CC). 

Recently, a QMC technique based on the resonating va- 
lence bond wavefunction was introduced in Ref . II 1 9fl . and fur- 
ther developed later. This approach represents a very efficient 
implementation of the valence bond Pauling idea, discovered 
by P.W. Anderson in the field of strongly correlated electrons 
112011 : the Jastrow Antisymmetrised Geminal Power (JAGP). 
This wavefunction has been demonstrated to be effective in 
describing highly correlated diatomic molecules like the C2 as 
well as 7r—7r interacting complexes il7lll9ll21ll . In the present 
article we present a Variational Monte Carlo (VMC) and Lat- 
tice Regularized Diffusion Monte Carlo (LRDMC) study of 
the water dimer dissociation energy and dispersion curve, us- 
ing as a variational ansatz the JAGP wavefunction. An im- 
portant advantage of the JAGP VMC approach resides in the 
possibility to dissect, in a simple way, the energy contribu- 
tions of the different terms composing the wave function, like 
dynamical electron correlations and the intermolecular cova- 
lent contribution. Dynamical electronic correlation associated 
to the charge fluctuations and van der Waals interactions, are 
indeed included in the wavefunction within by Jastrow terms, 
whereas static correlation are described by the resonance of 
valence bond singlets in the AGP. The amount of binding 
energy arising from the correlated dynamical charge fluctua- 
tions, related to the vdW forces, can be therefore estimated by 
the evaluation of the energy contributions of the Jastrow fac- 
tors. On the other side, following the Pauling idea of chemical 
bonding we can calculate the energy contribution of the inter- 
molecular covalent term and get insight on the covalent nature 
of the hydrogen bonding mechanism. 



II. COMPUTATIONAL METHODS 

Geometries. As nuclear coordinates of the water monomer 
we used the experimental equilibrium geometry I22I1 with 
an O-H bond length of 0.9572 A and a H-O-H angle of 
104.52°. For the dimer we used the linear configuration with 
C s symmetry, oxygen-oxygen distance 2.976 A I23I l24ll and 
0\ — H\ ■ ■ ■ O2 angle of 180°. We used the internal geome- 
try of each monomer as the experimental one. For the disper- 
sion curve we simply used the geometries obtained by shifting 
away the two monomers along the 0\ — Hi . ■ ■ ■ 0% binding 
axis and keeping fixed their relative orientation. Effects of nu- 
clear relaxation upon binding do not affect our estimations, as 
we verified by calculating the binding energy with the geom- 



etry from CCSD(T) calculations JH. 

Variational Monte Carlo and the JAGP wavefunction. 

As variational ansatz we use the JAGP wavefunction intro- 
duced in references 11911 and 12111 . The wavefunction ^jagp 
of a system of N electron is defined by the product of a sym- 
metrical Jastrow term J and an antisymmetrical determinantal 
part ^agp- 

^JAGp{n, —,Tn) = *AGp(ri, ...,rjv)J(ri,...,rjv) (1) 

The determinantal part ^>agp is the antisymmetrized product 
of spin singlets. The pairing function in singlet system with- 
out spin polarization is described by: 

*agp = A[*(rl,ri) • ■ • Hrl /2 ,r l N/2 )] (2) 

where A is the operator that antisymmetrizes the product of 
N/ 2 geminal singlets $(r T ,r+) = ip(r^ , r-L)L/\/2(| Tl> - 
I It))- The spatial part of the geminals are expanded over an 
atomic basis set: 

^(r\r l ) = ^ a , fc (rV) (3) 

a . b 

Vw^V) = XX'm<Mr T )<Km(r ; ) (4) 
1, m 

where the indexes I, m runs over different orbitals centered on 
nuclei a, b. The Jastrow factor J is further split into one-body, 
two-body and three-body terms (J = Ji 0/2^3)- The J\ and J2 
terms deal with electron-electron and electron-ion correlation 
respectively. The two body (one body) Jastrow depends only 
on the relative distance r^j = |rj — r,-| between each elec- 
tron pair (electron-ion pair) and has been parametrized 
by a simple function u(r.ij) = (1 — exp(— &ry))/2& that 
rapidly converges to a constant when ry became large. lfl7ll 
In this way the large distance behaviour of the Jastrow is de- 
termined only by the J3 Jastrow factor, that contains all vari- 
ational freedom left and in particular, as we shall see later 
on, the slowly decaying vdW correlations. Therefore we have 
chosen to parametrize this important part of our correlated 
wave-function in a systematic and exhaustive way, similarly 
to what we have done for the AGP contribution: 

J 3 (ri,...,rjv) = exp ^^$ J (r*, r j)j 

a, 6 

Kbi^j) = J2&aA^Xm(r 3 ) (6) 

Both the determinantal (j> a j and and Jastrow <j) J a l orbitals are 
expanded on gaussian basis sets centered on the correspond- 
ing nuclear centers a and b. By increasing the atomic basis 
set one can rapidly reach the "complete basis set limit" be- 
cause all cusp conditions are satisfied by an appropriate and 
simple choice of the J\ ( satisfying the electron-ion cusp) and 
J2 (satisfying the electron-electron cusp) terms. lfl7ll 
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All variational parameters, such as the Jastrow parameters, 
the {g} and the {A} matrices of equations [3] and [5] as well 
as the exponents and the coefficients of the gaussian orbitals 
have been optimized by e nerg y minimization following the 
methods described in refs. I261I27I1 . 

The oxygen valence-core interaction was described using 
the recently reported energy-consistent pseudopotentials J28ll . 
A VMC calculation for the dimer system with the larger ba- 
sis set and with 0.1 mH accuracy run for about twelve hours 
on eight AMD Opteron 280 CPUs at the CASPUR computer 
centre. Full wavefunction optimization was about a factor 4 
more time consuming. 

Diffusion Monte Carlo. A systematic way for improv- 
ing the quality of a variational wave function is to perform 
a Diffusion Monte Carlo calculation, filtering the ground state 
properties by a diffusion process B29I1 . Actually, due the pres- 
ence of the fermionic problem, the DMC is implemented 
within the fixed node (FN) approximation [30], by impos- 
ing that the final ground state has the same nodal structure 
of the trial WF. In this work we use a slightly modified ver- 
sion of the DMC method, the Lattice Reguralized Diffusion 
Monte Carlo (LRDMC). In this method, the continuum Monte 
Carlo moves are made by discrete finite steps defined by two 
lattice spaces a and a! . By using an incommensurate ratio 
a! I a the electronic trajectory fills the entire space, thus avoid- 
ing most lattice artifacts. The introduction of the lattice im- 
plies that there are a finite number of possible moves during 
the diffusion process, and this allows one to avoid the local- 
ity approximation and to restore the upper bound property of 
DMC.iUl [H] Within this regularization the exact Hamilto- 
nian H is replaced by a lattice reguralized one H a such that 
H a — > H for a — > 13111 . We used as lattice spaces the values 
a = 0.1, 0.2, 0.3, 0.5 a.u. and then the energy were extrapo- 
lated to zero lattice space. 

Since the dipole moment operator does not commute with 
the Hamiltonian, in LRDMC we evaluated the estimator fj, = 
2hlrdmC - VVMC, where ulrdmc is the LRDMC mixed 
average value extrapolated to a = 0. A LRDMC calculation 
for the dimer system with the larger basis set, 0.1 mH accu- 
racy and a = 0.2, run for about twelve hours on eight AMD 
Opteron 280 CPUs at the CASPUR computer centre. 

DFT Calculations. For the sake of comparison we perform 
DFT calculations using a plane wave basis set as implemented 
in the CPMD code 113 311 . For the exchange and the correlation 
part of the universal functional we used BLYP generalized 
grad ient corrections IL34T,l35ll and the hybrid functional B3LYP 
03611 . Core electrons were taken into account using norm- 
conserving Troullier-Martins type pseudopotentials yj . We 
also performed calculations with Dispersion-Corrected Atom- 
Centered Potentials (DCACP) [16] as described in reference 
113 811 . The Kohn-Sham orbitals were expanded in plane waves 
up to a cutoff of 125 Rydberg. 



VMC 


3-body Jastrow Basis 


Eh 2 


E (H 2 0) 2 


D e 


2s2p-local[0]ls[H] 


-17.2279(1) 


-34.4585(2) 


-0.0024(4)[- 1.5(0.3)] 


2s2p[0]ls[H] 


-17.2388(2) 


-34.4807(5) 


-0.0031(7)[-1.9(0.4)] 


2s4p[0]ls[H] 


-17.24089(5) 


-34.4874(1) 


-0.0056(2)[-3.5(0.1)] 


2s6p[0]ls[H] 


-17.24119(8) 


-34.4886(1) 


-0.0062(4)[-3.9(0.3)] 


2s6p[0]lslp[H] 


-17.2435(1) 


-34.4940(1) 


-0.0071(2)[-4.5(0.1)] 


LRDMC 


2s2p-local[0]ls[H] 


-17.2576(2) 


-34.5228(2) 


-0.0076(3) [-4.8(0.2)] 


2s2p[0]ls[H] 


-17.2613(1) 


-34.5303(1) 


-0.0077(2) [-4.8(0.2)] 


2s4p[0]ls[H] 


-17.2619(1) 


-34.5315(1) 


-0.0077(2) [-4.8(0.2)] 


2s6p[0]ls[H] 


-17.2619(1) 


-34.5314(1) 


-0.0076(2) [-4.8(0.1)] 


2s6p[0]lslp[H] 


-17.2620(1) 


-34.5318(1) 


-0.0078(2) [-4.9(0.1)] 



TABLE I: VMC and LRDMC energies for the water monomer 
and dimer (atomic units). The dissociation energy calculated 
as D e =E {H20)2 -2E H2 o is reported in the last column in atomic 
units and (in square brackets) in kcal/mol. 



III. RESULTS AND DISCUSSION 
A. Dissociation energy and charge fluctuations 

In this section we report our results on the water dimer at 
the experimental binding distance, and we investigate the in- 
fluence of different Jastrow terms of the wavefunction on the 
dissociation energy. In the pairing determinant the oxygen 
atoms are described using a gaussian basis set of 4s4p con- 
tracted to [ls2p], whereas we have only an uncontracted Is 
shell for the hydrogen. We verified that the inclusion of a d- 
wave shell does not affect the binding energy giving only a 
rigid shift of the total energy of the dimer and the monomer 
within LRDMC. 

On the contrary, more subtle effects have been observed 
in the structure of the three-body J3 Jastrow factor. This 
term includes in the wave function additional dynamical elec- 
tron correlations and contributes to the proper behavior of 
the electronic charge distribution. The correct description 
of the charge correlations reveals crucial for the inclusion of 
the vdW interactions, being originated by the correlations be- 
tween charge fluctuations in different spatial regions lf39ll . 

In table Q] we report the JAGP monomer energy, Eh 2 o, 
the dimer energy, Erjj 2 o) 2 > an d tne dissociation energy D e 
for increasing three-body Jastrow basis sets. The dissocia- 
tion energy of the water dimer has been calculated simply 
as D e = E dimer - 2E monomer . As the number of p-wave 
shells is increased, we observe an improvement of the bind- 
ing energy, eventually obtaining at VMC level a value of 
D e = 4.5(0. l)kcal/mol. 

The reported LRDMC results, extrapolated to the a = 
limit, appear to have a much faster convergence in term of to- 
tal and dissociation energies. This is due to the nature of the 
projection method which accuracy relies only on the nodal 
surface of the trial wave function. Our results indicate that the 
VMC optimized nodal surface, and therefore the correspond- 
ing DMC energy, is only slightly affected by the basis set ex- 
tension of the Jastrow factor. On the other hand, at the VMC 
level, the size of the Jastrow factor is crucial for improving the 
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VMC 


3B Jastrow Basis 


f/f 2 [Debye] 


M2[H 2 0] [Debye] 


2s2p-local[0]ls[H] 


2.116(17) 


2.805(20) 


2s2p[0]ls[H] 


1.935(12) 


2.834(23) 


2s6p[0]ls[H] 


1.880(8) 


2.692(14) 


2s6p[0]lslp[H] 


1.890(8) 


2.597(12) 


Extrapolated 


2s6p[0]ls[H] 


1.874(10) 


2.648(18) 


2s6p[0]lslp[H] 


1.870(10) 


2.603(13) 



TABLE II: Dipole moment of water monomer and dimer. VMC 
estimates are reported in the top part of the table. Extrapolated 
values are reported in the bottom part of the table. 



binding energy of the dimer. This is probably due to the role 
of the Jastrow in localizing charges and in introducing dynam- 
ical correlations. On this regard the role of p-wave shells in 
the binding energy will be discussed later in more details. 

The LRDMC results are obtained extrapolating at zero lat- 
tice space and give a binding energy of 4.9 ± 0.1 kcal/mol, 
in good agreement with the high level quantum chemistry cal- 
culations: Klopper et. al. have reported 4.99 kcal/mol and 
5.02(5) kcal/mol, as basis set limit for MP2 and Coupled 
Cluster calculations [25]. Recently a QMC study[18] has re- 
ported a value of 5.4(1) kcal/mol, obtained without optimiza- 
tion of the orbitals in the determinant, which is directly taken 
from a B3LYP calculation. 

Our LRDMC results are in the range of previously re- 
ported all-electron and pseudopotential QMC calculations 
ill 8[ l24l PtoTl . We also agree with experimental results al- 
though they suffer uncertainty due to theoretical estimation 
of the ZPE. Actually when compared with the experimental 
dimer dissociation energy D% xp , the difference between the 
zero point energy (ZPE) of the monomer and the dimer should 
be also taken into account: Df p = D - 2ZPE monomer + 
ZPE d%mer . The experimental energy reported hereafter is 
therefore corrected by this quantity calculated by theory or 
estimated by experiments [24]. As pointed out in Ref. lfl7ll 
the JAGP wavefunction is certainly size consistent for the two 
water monomers, only when the complete basis set limit is 
reached for the Jastrow factor. This property can be used to 
check the basis-set accuracy of the three-body Jastrow term. 
To verify the size consistency of the wave function we cal- 
culated the dissociation energy D* e by separating the two 
monomers at large distance. D* agrees with D e . 

In QMC calculations, correlation functions different from 
the energy are often very sensitive to the quality of a wave- 
function. We have therefore calculated the monomer and 
dimer dipole moment /1 that can be easily computed at the 
VMC level and at LRDMC level using the mixed estimator. 
In the case of the water monomer both the variational and 
the LRDMC dipole moments are rather close to the experi- 
mental value of 1.8551? [3]. The correction introduced by 
LRDMC is a slight downshift of the extrapolated estimator, 
/x = 1.870(10) D, in agreement with other QMC calculations 
lfl8h and ab initio methods ll4lll . 

We turn now the attention to the role of dynamical correla- 



tions included through the Jastrow term. As discussed above 
the inclusion of p-wave orbitals in the J3 Jastrow term has 
a significant effect on the binding energy. Similarly, the ef- 
fect of the J3 basis set is also visible on the dipole moments 
in table [II] One reason of this influence can be attributed to 
electrostatic interactions, since the 3-body term is important 
for the charge distribution. Another relevant effect of the J3 
term is the modulation of the van der Waals interactions. Van 
der Waals forces are indeed quantum mechanical effects orig- 
inating by the interaction between instantaneous dipoles, or, 
using a second order perturbation theory perspective l42tl . by 
the correlated transition of a couple of electrons from occu- 
pied to unoccupied states. Given thus two atomic centers a 
and b at large distance the J3 term in Eq. [5] can be expanded 
for small value of gfj and then applied to a single gemi- 
nal pair (see Eqf3), ^a,b(T,r^). The result can be viewed 
as a correlated transition of two electrons located in different 
atomic centers from occupied orbitals to unoccupied orbitals 
with higher angular momentum[42[. More generally the ef- 
fect of the J3 at large distance has the same structure of the 
vdW perturbative term if on each atomic center the basis used 
for the Jastrow contains odd orbitals with respect to the spatial 
reflection, namely when the Jastrow basis set contains at least 
p wave orbitals. In principle a small vdW contribution can 
derive also from high angular momentum orbitals included in 
the geminal expansion. In this work however, in order to dis- 
entangle the genuine vdW contribution, we have avoided to 
use polarization orbitals in the AGP, that, as discussed before, 
do not affect the binding energy. In this way the instantaneous 
correlated polarization induced by the J3 term allows to in- 
clude vdW interactions in a transparent variational form. 

To understand the effects of the J3 terms on the dissociation 
energy, we calculated the variational energy of the wavefunc- 
tion obtained excluding inter-molecular g?' terms in Eq. [5]as 
reported in table [III] In particular we considered the H — O 
and — contributions in the p-p channel and eventually we 
eliminated all intermolecular terms (last row of the table [Tllb. 
Data show non-additivity of the energy loss, as expected by 
interactions arising from polarization effects [43]. Among the 
p — p wave contribution, the oxygen-oxygen channel seems to 
be the most relevant term in the Jastrow expansion. It is worth 
noting that the total dipole moment of the dimer depends only 
weakly on the intermolecular J3 Jastrow terms, see table [III] 
This indicates that the distribution of the electronic charge is 
not much affected by the missing terms. The energy differ- 
ences are then due to the part of the dynamical correlation 
involving correlated excitations to p states. The energy loss 
in the binding energy can therefore be attributed within our 
formalism to van der Waals interactions. 

It is of interest to compare our result to previous cal- 
culations based on symmetry-adapted-perturbation-theory 
(S APT) H44l I45H that estimated the contribution of dispersion 
forces to the water dimer hydrogen bond. This contribution 
amounts to about -1.75 kcal/mol as reported in table V of Ref. 
j45ll . Albeit the energy is not partitioned the same way in 
the two approaches, the assessment given by SAPT in good 
agreement with our estimation of -1.5(2) kcal/mol. 
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Pairing terms in J3 


E(H 2 o) 2 (a-u) 


AE (kcal/mol) 


Hh 2 o) 2 [D] 


Ful1 i9l,m} matrix 


-34.4940(1) 


0.0 


2.597(12) 


(Plnm (PlnJ)2 




+0.2(1) 


T AT 1 / 1 Ti 
L.KiLYyYL) 


(p[0])i (p[H]) 2 


-34.4938(1) 


+0.1(1) 


2.623(12) 


(p[H])j (p[0]) 2 


-34.4935(1) 


+0.3(1) 


2.610(13) 


(p[0])i (p[0]) 2 


-34.4918(1) 


+1.4(1) 


2.628(12) 


intermolecular p-p gf'^ =0 


-34.4916(3) 


+1.5(2) 


2.637(13) 



TABLE III: VMC energy and dipole moment of the water dimer 
for different Jastrow J 3 terms. The different J 3 are obained 
by canceling the p - p electronic correlation between atomic 
centers belonging to different molecules. The atomic center of 
the p wave is indicated between square brackets and the water 
molecule index is indicated by the pedex (1 or 2). The energy 
difference A.E with respect the complete g"^ matrix, first line, 
is also reported in kcal/mol. 
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B. Dispersion curve 



The VMC and LRDMC dispersion curve of the water dimer 
is reported in FigQ~|\. It has been calculated by computing the 
total dimer energy as a function of the oxygen-oxygen dis- 
tance without changing the internal geometry and the relative 
orientation of the monomers. We used the 2s6p[0] ls[H] basis 
set for the 3-body Jastrow, which, as reported before, guaran- 
tees size consistency during the dissociation process at large 
distances. 

The attractive tail of the water-water interaction potential 
is dominated by a dipole-dipole interaction energy. A poly- 



2{H 2 0) 



d a with 



nomial fit for d > 3.5 A shows that E : 
a = 3.2 - 3.3 for the VMC and the extrapolated LRDMC 
curves, respectively. 

The behavior of the dispersion curve at short distances is 
shown in the inset of Fig. [TJ\, together with a fit performed 
using a Morse potential. At the VMC level, the minimum 
of the curve as obtained by the fitting procedure, is at dis- 
tance d = 3.037(4) A, which is slightly shifted with re- 
spect to d = 2.976 A reported by experiments [23]. How- 
ever, it should be noted that, considering the error bars, the 
curve results to be rather flat around the equilibrium distance. 
LRDMC with a lattice space a = 0.2 a.u., and the LRDMC 
extrapolation to zero lattice space a — > 0, improve the loca- 
tion of the equilibrium distance. In this latter case the fitted 
minimum is at d = 2.982(1) A, which is very close to the 
experimental value. In figureQJi we report a comparison with 
pure or empirically parametrized Density Functional methods 
and symmetry-adapted perturbation theory (SAPT)|44, 46]. 
Data show that pure BLYP and B3LYP curves underestimate 
the dissociation energy whereas calculations performed with 
empirically parametrized DCACP pseudopotentials are closer 
to our LRDMC curve. SAPT curve is on the top of our results 
with small discrepancies at very short distance. 
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FIG. 1: (Color on line)Water dimer dissociation. The total en- 
ergy of the water dimer is reported as a function of the oxygen- 
oxygen distance. In Panel A the VMC and LRDMC results are 
reported. In the inset graph the behavior around the minimum 
is zoomed in. In Panel B the LRDMC curve is compared with 
other methods. SAPT values have been taken from ref. [46]. 



C. Covalent contribution to hydrogen bonding 

In the proximity to the equilibrium distance the interplay 
between electrostatic and pure quantum effects is expected to 
be relevant. Although a unique and commonly accepted defi- 
nition of covalent contribution in hydrogen bond is still miss- 
ing, within the formalism of the JAGP wavefunction, we can 
define the covalent energy contribution as the energy contribu- 
tion arising from the intermolecular pairing terms of the RVB 
determinantal part of the wavefunction (see Eq. |3J- 

The "chemical bond" between the two molecules is indeed 
due to the superposition of all singlet terms in the geminal 
expansion that connect two nuclei belonging to different wa- 
ter molecules. This is schematically illustrated in the panel 
A of Fig. |2 In order to evaluate the covalent contribution 
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FIG. 2: (Color on line)Top panel: Pictorial view of the inter- 
molecular pairing term of the determinant part of the wave- 
function. In Panel A all the intermolecular pairing are drawn. 
In Panel B we represent the wavefunction with the intermolec- 
ular pairing term set to zero. In the bottom part of the Figure 
the total energy of the water dimer is reported as a function of 
the oxygen-oxygen distance between the two monomers. Two 
wavefunctions are compared, one with all the intermolecular 
pairing term optimized (black circle) and the other one with in- 
termolecular pairing terms set to zero (red square). 



we proceed as follows. We cut the intermolecular pairing va- 
lence bonds in the pairing function, by imposing ^ 0) & = if 
a and b belong to different monomers (as sketched in panel 
B of Fig|2j. Then the wave function is re-optimized with the 
above constraint in order to correctly include the electrostatic 
effects and the slowly decaying vdW correlations present in 
our Jastrow factor. In Fig. [2] C we report, as a function of 
the oxygen-oxygen distance, the binding energy calculated 
with the full wavefunction (circle) and with the wave function 
lacking the intermolecular valence bond terms (square). The 
difference between the two curves vanishes as the molecules 
reach a 0-0 distance of 3.5 A. We point out that, by cutting 
the intermolecular pairing terms, the minimum of the energy 
dispersion slightly shifts to a larger equilibrium distance. 

At the equilibrium distance we found that the contribution 
of the intermolecular pairing terms, computed at the VMC 
level, is Ai nter = 1. 1(0. l)kcal/mol, corresponding to about 
24% of the computed dimer binding energy. Our estimate of 



the covalent contribution, defined above, can be compared to 
what other authors found using different theoretical frame- 
works and different definition and that is generally referred 
to as intermolecular charge transfer (CT) l)47l lil l49l l50ll. 
In the seminal works based on Morokuma decomposition of 
the binding energy ll47l l48ll the CT contribution to hydrogen 
bonding is estimated in the range — 1.3 -i — 1.8 kcal/mol., thus 
about 25% of the total binding energy. A slight smaller contri- 
bution, about 11% resort from block-localized wave function 
approach proposed by Mo et al. 14911 . They also reported that 
CT contribution vanishes at about 3.5 A in good agreement 
with our finding. It is interesting to observe that the damping 
of the oscillation at d ~ 4 A of the Fourier Transforms of the 
Compton profile, has been interpreted in ref. [ 12] as a cut-off 
for the covalent contribution. ||9l|l2|]. However such an inter- 
pretation of the experiments is not fully accepted [ 1 Q UI ill . 



IV. CONCLUSIONS 

The understanding of hydrogen bond systems is still a chal- 
lenge for computational chemistry. Even for small molecu- 
lar systems the weakness of the interactions and the critical 
role of electron correlation require to use affordable correlated 
quantum chemistry methods. The interplay between interac- 
tions different in nature, such as dispersion forces and inter- 
molecular charge transfer are in many cases crucial for the 
proper description of the bond properties. 

We have shown that the Quantum Monte Carlo method is 
effective for describing the hydrogen bond between two water 
molecules. The calculated binding energy matches the exper- 
imental value and the estimates from other advanced method- 
ologies. Good agreement with experiments is also achieved 
for the computed dipole moments. Thanks to the good size 
scaling properties and the embarrassingly parallelism of QMC 
algorithms, these methods appear extremely competitive in 
the context of massive parallel computation. 

Moreover, some conceptual advantages rely on the struc- 
ture of the AGP wavefunction, a correlated valence bond rep- 
resentation of the electronic system. The AGP formalism 
gives the possibility to work back on an intuitive picture of 
localized chemical bonds such as the Pauling's superposition 
of Lewis structures. Thanks to the fully correlated structure 
of the wavefunction, this picture can be used without compro- 
mises in term of accuracy. 

Upon interpretation of the wavefunction terms, we estimate 
at the VMC level the covalent contribution to account for 
1.1(2) kcal/mol. A similar contribution to the binding energy 
is given by correlated dipolar vdW fluctuations that account 
for 1.5(2) kcal/mol. 

The quality of our results on the water dimer encourages the 
application of the method to larger hydrogen bonded systems 
such as water clusters or small biomolecules. To reduce the 
computational costs it would be desirable to keep down the 
number of variational parameters when the size of the system 
increase. In this respect, different strategies are under investi- 
gation. 
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